#include "config.h"

#include "lstatfit.h"

#define GAP	0.00001

ecm_double
lognormal_goodness(msl_fit_t *fit)
{
	ecm_uint32	i;
	ecm_double	*xvalues = fit->xvalues;
	ecm_double	*yvalues = fit->yvalues;
	ecm_double	mu, sigma;
	ecm_double	sq_diffs = 0;

	mu = gsl_vector_get(fit->s->x, 0);
	sigma = gsl_vector_get(fit->s->x, 1);

	for (i = 0; i < fit->n; i++) {
		double	t = (xvalues[i] + GAP) / (100 + 2 * GAP) * 100;
		double	front_part, back_part;
		ecm_double	prob;

		front_part = 1 / (t * sigma * sqrt(2 * M_PI));
		back_part = exp(-(pow(log(t) - mu, 2)) / (2 * sigma * sigma));

		prob = front_part * back_part;

		sq_diffs += pow(yvalues[i] - prob, 2);
	}

	return sq_diffs;
}

void
lognormal_dist_report(FILE *fp, msl_fit_t *fit, ecm_uint32 xdiff)
{
	ecm_uint32	i;
	ecm_double	*xvalues = fit->xvalues;
	ecm_double	mu, sigma;

	mu = gsl_vector_get(fit->s->x, 0);
	sigma = gsl_vector_get(fit->s->x, 1);

	ecm_fprintf(fp, "# mu:%lf sigma:%lf\n", mu, sigma);
	
	for (i = 0; i < fit->n; i++) {
		double	t = (xvalues[i] + GAP) / (100 + 2 * GAP) * 100;
		double	front_part, back_part;
		ecm_double	prob;

		front_part = 1 / (t * sigma * sqrt(2 * M_PI));
		back_part = exp(-(pow(log(t) - mu, 2)) / (2 * sigma * sigma));

		prob = front_part * back_part;

		ecm_fprintf(fp, "%lf %8.7f\n", t + xdiff, prob);
	}
}

int
lognormal_dist_f(const gsl_vector *x, void *data, gsl_vector *f)
{
	msl_fit_t	*fit = (msl_fit_t *)data;
	ecm_uint32	n = fit->n;
	ecm_double	*xvalues = fit->xvalues;
	ecm_double	*yvalues = fit->yvalues;
     
	double	mu = gsl_vector_get(x, 0);
	double	sigma = gsl_vector_get(x, 1);

	ecm_uint32	i;
     
	for (i = 0; i < n; i++)	{
		double	t = (xvalues[i] + GAP) / (100 + 2 * GAP) * 100;
		double	front_part, back_part;

		front_part = 1 / (t * sigma * sqrt(2 * M_PI));
		back_part = exp(-(pow(log(t) - mu, 2)) / (2 * sigma * sigma));

		double	Yi = front_part * back_part;
		gsl_vector_set(f, i, Yi - yvalues[i]);
	}

	return GSL_SUCCESS;
}

int
lognormal_dist_df(const gsl_vector *x, void *data, gsl_matrix *J)
{
	msl_fit_t	*fit = (msl_fit_t *)data;
	ecm_uint32	n = fit->n;
	ecm_double	*xvalues = fit->xvalues;

	double	mu = gsl_vector_get(x, 0);
	double	sigma = gsl_vector_get(x, 1);

	ecm_uint32	i;
     
	for (i = 0; i < n; i++)	{
		double	t = (xvalues[i] + GAP) / (100 + 2 * GAP) * 100;
		double	dmu, dsigma;
		double	dmu_part;
		double	ds_front_part, ds_back_part;
		double	front_part, back_part;

		front_part = 1 / (t * sigma * sqrt(2 * M_PI));
		back_part = exp(-(pow(log(t) - mu, 2)) / (2 * sigma * sigma));

		dmu_part = (log(t) - mu) / (sigma * sigma);
		ds_front_part = -1 / (t * sqrt(2 * M_PI) * sigma * sigma);
		ds_back_part = back_part * (pow(log(t) - mu, 2) / pow(sigma, 3));

		dmu = front_part * back_part * dmu_part;
		dsigma = ds_front_part * back_part + front_part * ds_back_part;

		gsl_matrix_set(J, i, 0, dmu);
		gsl_matrix_set(J, i, 1, dsigma);
	}
	return GSL_SUCCESS;
}
